% simulation of ddot(x)=-k sign(x)
params=struct();
tfinal=40;   
params.k=@(t)1;
initial=[2;2];
[d0,s0]=rhs(0,initial,params);
profile=repmat(params,4,1);
profile(1).k=@(t)1;
profile(1).step=0.01;
profile(1).title="const1";

profile(2).k=@(t)exp(-t);
profile(2).step=0.01;
profile(2).title="exp-1";

profile(3).k=@(t)abs(sin(t));
profile(3).step=0.1;
profile(3).title="abssint";

for ip=1
    params=profile(ip);
    [t,x]=rk4(@(t,x)rhs(t,x,params),0,initial,200,params.step);
    s=repmat(s0,length(t),1);
    for i=1:length(t)
        [di,si]=rhs(t(i),x(i,:)',params);
        s(i)=si;
    end

    % plot(t,x(:,1))
    % plot(t,x(:,2))
    tiledlayout(3,1)
    nexttile;hold on;
    semilogy(t,abs(x(:,1)))
    ylabel("x")
    nexttile;hold on;
    semilogy(t,abs(x(:,2)))
    ylabel("v")
    nexttile;hold on
    plot(t,horzcat(s.u))
    legend(["u"])
    exportgraphics(gcf,"out/"+params.title+".png")
    
    figure();
    plot(x(:,1),x(:,2))
    xlabel("x")
    ylabel("v")
    axis equal
    grid on
    exportgraphics(gcf,"out/"+params.title+"traj.png")
end

function [dxdt,s]=rhs(t,states,p)
    s=struct();
    x=states(1);
    v=states(2);
    s.u=-p.k(t)*sign(x);
    dxdt=[v;s.u];
end